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Abstract 

Parallel tempering (PT) Monte Carlo simulations have been applied to a variety of systems 
presenting rugged free-energy landscapes. Despite this, its efficiency depends strongly on the tem- 
perature set. With this query in mind, we present a comparative study among different tempera- 
ture selection schemes in three lattice-gas models. We focus our attention in the constant entropy 
method (CEM), proposed by Sabo et al. In the CEM, the temperature is chosen by the fixed 
difference of entropy between adjacent replicas. We consider a method to determine the entropy 
which avoids numerical integrations of the specific heat and other thermodynamic quantities. Dif- 
ferent analyses for first- and second-order phase transitions have been undertaken, revealing that 
the CEM may be an useful criterion for selecting the temperatures in the parallel tempering. 
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I. INTRODUCTION 

Enhanced sampling tempering approaches, such as the parallel tempering (PT) [l| and 
simulated tempering (ST) [2j, have played an important role in the study of complex systems 
such as spin glasses |3j|, protein folding |4|, biomolecules and others. The basic idea consists 
of using the "information" obtained at high temperatures to systems at low temperatures, 
allowing the system to escape from metastable states and providing an appropriate visit of 
the configuration space. 

Despite this, its efficiency depends strongly on the temperature set. In the last years, 
different procedures have been proposed for the choice of both the temperature set and the 
number of replicas. Kofke |5j has related the average acceptance probability of replicas p acc 
with its difference of entropy. Predescu et al. [6], based on the link between acceptance rate 
and specific heat, argued that the optimal distribution of intermediate temperatures should 
follow a geometric progression if the specific heat is nearly constant. Kone et al. [7] proposed 
that temperatures in the parallel tempering should be chosen such a way that about 20% 
of exchanges are accepted, when the specific heat is also constant. Katzgraber et al. 8] 
considered a feedback optimized parallel tempering Monte Carlo that gives the temperature 
set by measuring the replica diffusion through the temperature space. More recently Sabo 
et al. [9] introduced the constant entropy method (CEM), where the temperature of replicas 
is chosen provided the difference of entropy is held fixed. An advantage of both CEM 
and feedback optimized methods is that near the criticality temperatures become more 
concentrated, what increases the frequency of the replica exchanges when compared with 



other schemes. Results for the Ising model 



9] revealed that they are equally efficient and 



superior than other schemes. However, the CEM seems to be computationally simpler than 
the method by Katzgraber et al. 

In this paper, we give a further step in order to ascertain the role of the temperature 

schedule in the parallel tempering. We consider a comparative study among five different 

schemes for selecting the temperatures to be used: arithmetic and geometric distribution 

I I 

of temperatures, arithmetic distribution in inverse temperature [U[ , ad — hoc distribution 
of temperatures and the CEM. In the ad hoc ensemble, temperatures are chosen in such a 
way that exchanges between adjacent replicas are about 30% 0, n|. In order to select the 



temperatures according to the CEM, we use the transfer matrix method 



121 ]. in which the 



entropy is obtained directly from Monte Carlo simulations, hence avoiding numerical inte- 
grations of thermodynamic quantities, such as the specific heat and energy. We shall apply 
;he comparative study to different lattice models undergoing phase transitions such as Ising 
13 1, Blume- Emery-Griffiths (BEG) Q and Bell-Lavis (BL) water model JjJ. Different 



analyses have been undertaken and for all schemes considered here the CEM revealed to be 
more advantageous. 

This paper is organized as follows: In Sec. II we describe the approaches employed here, 
in Sec. Ill we present the models, in Sec. IV we discuss the numerical results and in Sec. V 
we give the conclusions. 

II. PARALLEL TEMPERING AND THE CEM 

The basic idea of the PT is that configurations from high temperatures are used to 
perform an ergodic walk in low temperatures. To this end, one simulates simultaneously a 
set of R replicas ranged in the interval {T 1; ...,Tr} by means of a Metropolis like algorithm 



161 ] (although in general the interest lies in the lowest temperature T = T\). The actual MC 
simulation is composed of two parts: In the first part, a given site k of the replica % is chosen 
randomly and its variable o~k may change to a new value a' k according to the Metropolis 
prescription p { = min{l, exp[— ^AH]}, where AH = H(cr') — H{(t) is the energy change 
due to the transition and j3i = (/c^Ti) -1 . In the second part, arbitrary pairs of replicas (say, 
at Tj and 7} with microscopy configurations a' and a") can undergo temperature switchings 
with probability 

Vi^j = min{l, exp[(A - P j ){H{a f ) - U{o"))\}. (1) 

In principle, the replicas i and j may be adjacent or not. Usually one adopts only ad- 
jacent exchanges, since the probability of a given swap decreases by raising T. In some 



cases, however, non-ad t 
from metastable states 



acent exchanges have revealed essential for the system to escape 
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2l|. 



As mentioned in the Introduction, we are going to consider different procedures for deter- 
mining the temperature set. The former schemes (both arithmetic and geometric schedules) 
are rather simple, since they require only the extreme temperatures for determining the 
whole set. In the ad hoc distribution, one starts from T\ for a given R and determines the 
R—l temperatures in such a way that the exchange probability between adjacent replicas is 



about 30%. The CEM [9] consists of adding intermediate temperatures with fixed difference 
of entropy. More specifically, given the extreme temperatures 7\ and Tr with entropies per 
volume Si and s#, respectively, we add R — 2 intermediate temperatures Tj whose entropy 
Si is Sj = si + (i — 1) x As, where As = (sr — Si)/R — 1. Each value of s will be calculated 
through the thermodynamic equation s = where u = (H) and / is given by 

The transfer matrix method, used for obtaining / (and consequently s), is implemented by 
dividing a lattice with V sites in N layers with L "spins" (V = L x N). The associated 
Hamiltonian is given by 

H = Y,H(S k ,S k+1 ), (3) 

k=l 

where S k = (o- ltk ,a 2 ,k, ■■■, cr L,k), and due to the periodic boundary conditions SV+i = Si- 
The probability P(S\, S 2 , Sn) of a given configuration is given by 

P{S 1 , S 2 , S N ) = ^T(Si, S 2 )T(S 2 , S 3 )...T(S N , S 1 ), (4) 

where T(Sk, Sk+i) = exp(— /31-L(Sk, Sfc+i)) i s an element of the transfer matrix T and Z = 
Tr( T^) is the partition function. The marginal probability distributions P(Si) and P(Si, S 2 ) 
are given by 

P(S 1 ) = ±T N (S U S 1 ), (5) 



and 



^(Si, S 2 ) = ^T(S h S 2 )T N -\S 2 , Si). (6) 



We can use the spectral development of the matrix T given by 

T(S u S 2 ) = J2MSi)Xk<Pl(S 2 ), (7) 

N 

where 4>k{Si) is the normalized eigenvector of T and is the corresponding eigenvalue, to 
write the expressions 

^ = £Af, (8) 
k 

p(s 1 ) = WmSi)\UI(s 1 ), (9) 
L k 



and 



p(s 1 ,s 2 ) = ^T(s 1 ,s 2 )J2MS2)\t 1 <l>l(s 1 ). 

z k 



(10) 



In the limit N -> oo, Eqs. © 



and (HDJ become Ag , P(Si) = fc (5'i)^( i S 1 ) and 



P(Si, S 2 ) = j-T(Si, S 2 )(j)k(S 2 )(j)l(Si), respectively. Putting S2 = in the above expression, 



we arrive at 



A 



(T(S 1 ,S 1 )) 



(11) 



(fei,5 2 ) 

where A is the largest eigenvalue of T and / becomes / = — -^-lnA . The quantity Ss lt s 3 
is the Kronecker delta for Si and 5*2 which is equal to 1 when layers Si and S 2 are equal 
and zero otherwise. In the next section, we will write down the transfer matrix T(Si,Si) 
for the models studied here. Since the averages described above are evaluated over all N 



layers, from now on we are going to replace T 



transfer matrix method are found in Refs. 



[SjuSj) by T(Sk, Sk). More details about the 
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III. MODELS 



A. Ising model 

The Ising model [l3j is defined as follows: each site of the lattice is attached by a spin 
variable o that takes the values (7j = ±1 according to whether the spin is "up" or "down" 
(or equivalently, an occupied or empty site in a fluid jargon). The Hamiltonian is given by 



<i,j> i 



(12) 



where J is the interaction energy between two nearest neighbor spins and H is the magnetic 
field. For H = and low temperatures, the system presents a phase coexistence between 
two ferromagnetic phases which ends at the critical point T c = 2.269. ..,where T = ksT/J. 
The transfer matrix diagonal elements, which will be used for obtaining the entropy in the 
CEM, reads 



T(S k , S k ) = exp 



(3 { -J (1 + 0i,k <7i+i,k) + H a i)k ] 



i=l 



(13) 



where the sum is performed over a layer with L sites. 
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B. BEG model 

The BEG model [14J] is a generalization of the Ising model, where each site is allowed to 
be empty or occupied by two distinct species. It is defined by the Hamiltonian 

n = - £ iM *i + Ka 2 ^) + °l (14) 

<i,j> i 

where Oi = 0, if the site i is empty and ±1 if i is occupied by one of the species. Parameters 
J and K are interaction energies and D denotes the chemical potential. The BEG model 
displays a rather rich phase diagram, whose features depend on the ratio K = Kj J. As far as 
the regime K > —0.5 and low temperatures is concerned, the system displays liquid (p 7^ 0) 
and gas phases (p = 0) for high and low chemical potentials, respectively (or equivalently for 
low and high values of D = D/J, respectively). For T — 0, the liquid-gas phase coexistence 
takes place at D* = z(K + l)/2, where z is the coordination number. The transfer matrix, 
that will be used for evaluating the entropy, is given by 

L 

T(S k , S k ) = exp [(3^2 ( J a i+i,k&i,k 

i=l 

+(J-D + K(l + a 2 +l!k ))al k )]. (15) 

C. Bell-Lavis model 

The Bell-Lavis (BL) model is defined on a triangular lattice where each site may be empty 
(<x; = 0) or occupied by a water molecule (cij = 1). Each particle has two orientational states, 
that may be described in terms of bonding and inert "arms" rf , which take the values t\ 3 = 
or tI° = 1 when the arm is inert or bonding, respectively. Two nearest neighbor molecules 
interact via van der Waals e v d w and hydrogen bond energies whenever they are adjacent 
and point out their arms to each other {t^t? 1 = 1), respectively. The BL model is defined 
by the following Hamiltonian 

H = ~ a i °i ( e hb Ti if + e vdw ) -/J^tTj, (16) 

<i,j> i 

where p is the chemical potential. The BL model also displays a rich phase diagram, whose 
features depend on the ratio ( = e v d w /^hb- In particular, for ( = 0.1 one has three phases, 



denoted gas, low-density-liquid (LDL) and high-density-liquid (HDL) 



15 



As in the 



BEG model, the gas phase is devoided by molecules, whereas the LDL phase is characterized 
by a honeycomb like structure, with density of particles p and hydrogen bonds per molecule 
Phb given by p = 2/3 and phb = 3/2, respectively. In the HDL phase, the lattice is filled by 
molecules, p — 1, and the hydrogen bonds density per molecule is also given by p^ = 1- For 
low values of p = p/ehb, the system is constrained in the gas phase. By increasing p for fixed 
low T, a first transition between the gas and the LDL phase occurs. By increasing further p 
a second transition, from the phase LDL to the HDL, takes place. At T — 0, both transitions 
are first-order and occurs at p* = —3 (1 + C)/2 and p* = —G(, respectively. For T ^ the 



24]. 



former phase transition remains first-order, whereas the latter becomes second-order 
For ( = 0.1, the second-order and first-order lines meet in a tricritical point. 

The transfer matrix that will be used for evaluating the entropy in the CEM is given by 

L 

T(S k} S k ) = exp [XX°"M: (°"*>* + 2<J i+i,k) (17) 

i=l 

{tvdw + £hb T~i,k T i+l,k + A*)} 



IV. NUMERICAL RESULTS 



A. Ising model 



We have performed numerical simulations of the Ising model in a square lattice of size 
L = N = 20 using periodic boundary conditions. We have discarded 3 x 10 5 MC steps, 
in order to equilibrate the system and 3 x 10 6 MC to evaluate the appropriate quantities. 
In Fig. HJ we compare the entropy s/ks obtained via the transfer matrix method with 



exact results 22] . We have considered a set of R = 21 replicas ranged from T\ = 0.1000 to 
T21 = 10.0000 distributed according to the CEM. The agreement between results support 
the adequacy of the present procedure for obtaining s, which will be used for estimating the 
temperature set. 

In Table H] we compare our temperature estimates with those obtained by Sabo et al. {^J 
by means of numerical integration of the specific heat C v . The data are shown to be in 
good overall agreement, even though some small discrepancies can be observed at certain 
temperatures. They may be explained by either numerical uncertainties in the transfer 
matrix method, or by uncertainties in the integrations of C v , or by derivations of u in 



l 1 i 1 r 
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FIG. 1. Entropy per site s/ks versus temperature T for the Ising model. We have used L = N = 20 
and CE intervals. The result obtained via the transfer matrix method (stars) is compared with 
the exact result (solid line). 

the exact method, or even by all these sources together. In the first comparison among 
temperature schedules, we investigate the decay of time-correlation displaced functions C q 
at the critical point. This study is motivated by the fact that numerical simulations of 
second-order phase transitions via conventional algorithms are affected by a slow decay of 
C q (critical slowing down). On the other hand, cluster algorithms 25j reduce drastically 
this effect. This suggests that the analysis of C q may be a good measure for the comparison 
of different criteria used in the PT. The auto-correlation function C q of a given quantity q 
at the time r is given by 



C q (r) = ((q(t)-q)(q(t + r)-q))/a 



(18) 



where q is the mean value and a q the variance, respectively. In Fig. [2] we plot C q for the 
thermodynamic quantities m = (<Ji) (magnetization per site) and u = (7i) (total energy 
per site) for all distributions described above. We have considered R = 6 replicas ranged 
from T\ = 2.269 to T4 = 3.660, whose temperature set is showed in Table HT1 for the CEM 
and ad hoc distributions. By putting the temperature set into an array, we have considered 
exchanges between every adjacent and non-adjacent (here between every second and third) 
temperatures. 

Although all schemes give equivalent results for the steady state (for u and m), the 
quantity C q decays faster within the CEM. As it will be shown later, a similar behavior is 
verified when one calculates C q at the critical point for the BL model. In particular, by 
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CEM (exact) 


CEM-S (exact) 


CEM (20) 


CEM-S (20) 


0.1000 


0.1000 


0.1000 


0.1000 


1.4683 


1.4688 


1.4551 


1.4635 


1.6867 


1.6866 


1.6769 


1.6820 


1.8360 


1.8373 


1.8289 


1.8332 


1.9524 


1.9538 


1.9457 


1.9496 


2.0464 


2.0478 


2.0379 


2.0433 


2.1236 


2.1250 


2.1107 


2.1201 


2.1865 


2.1879 


2.1696 


2.1836 


2.2363 


2.2374 


2.2200 


2.2386 


2.2697 


2.2702 


2.2681 


2.2883 


2.3048 


2.3064 


2.3170 


2.3365 


2.3580 


2.3603 


2.3711 


2.3896 


2.4288 


2.4321 


2.4374 


2.4518 


2.5208 


2.5265 


2.5242 


2.5341 


2.6408 


2.6473 


2.6401 


2.6464 


2.8007 


2.8115 


2.7977 


2.8019 


3.0219 


3.0372 


3.0170 


3.0220 


3.3481 


3.3708 


3.3393 


3.3483 


3.8852 


3.9414 


3.8687 


3.8847 


4.9930 


5.1247 


4.9505 


4.9906 


10.000 


10.000 


10.000 


10.000 



TABLE I. Temperature CEM set for the Ising model obtained from the present approach (first 
and third columns) and by Sabo et. al [9J] (second and fourth columns). 



allowing only adjacent exchanges, both C m and C u also decay faster with the CEM than 
other schedules. 
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CEM ad hoc 
2.269 2.269 
2.345 2.400 
2.464 2.580 
2.658 2.830 
2.996 3.170 
3.660 3.660 

TABLE II. Temperature set of R = 6 replicas for the Ising model obtained for the CEM and ad 
hoc distributions. 




FIG. 2. Time displaced auto-correlation functions C m and C u versus time r (in MC steps) for the 
Ising model at the critical point for N = L = 20 and different schemes of temperature selection. 



B. BEG model 



We have performed numerical simulations for the BEG model in a square lattice of size 
L = N = 20 using periodic boundary conditions. We focus on the analysis for K = 3, H = 
and the low temperature T\ = 0.5000. In this case, a first-order phase transition between 



the liquids and the gas phase takes place at D* = 8.0000(1) 19] . Since the probability 
distribution at discontinuous transitions exhibit two peaks (corresponding to each phase), 
conventional Monte Carlo algorithms are not efficient at low temperatures, since the system 
requires a long time to pass from one peak to the other. In extreme cases, the peaks are 
separated by very high barriers and the system may get trapped in a given phase along the 
whole simulation and in this case it will be not ergodic. With these concepts in mind, we 

10 



0.9 
0.8 

P 

0.7 
0.6 

°' 5 2 4 6 8 10 

t(10 4 MC steps) 

FIG. 3. Time decay of the order parameter p versus time (in MC steps) for the BEG model at 
the phase coexistence point (D*, T\)= (8. 0000, 0.5000) for N = L = 20 and different distribution of 
temperatures. The tie line for p = 2/3 denotes its stationary value po- 

consider two analyses at the phase coexistence: the time evolution of the order parameter 
P = (°f ) toward its equilibrium value po = 2/3 |29j starting from a non typical configuration 
and the tunneling between the phases at the coexistence after discarding sufficient MC steps. 
The latter study will be carried out by measuring the fluctuation of p around po, since the 
trapping of the system in a given phase or in a mestastable state is expected to be signed 
by no relevant change of p. We have distributed temperatures between 7\ = 0.5000 and 
T R = 2.0000, whose entropies per site are given by Si/ks = 5 x 10~ 6 and sn/k B = 0.4971(1), 
respectively. By using in all set of R — 6 replicas, the CEM criterion leads to the 

intermediate temperatures T 2 = 1.5550, T 3 = 1.7650, T4 = 1.8780 and T 5 = 1.9400. As for 
the Ising model, we have considered exchanges between every adjacent and non- adjacent 
(here between every second and third) temperatures. In Fig. [31 we plot the time evolution 
of p starting from a lattice filled with particles. Note that by choosing the temperatures 
according to the CEM, the convergence of p toward po is faster than for other schemes. 
Although the results obtained for the ad hoc and CE cases are close, only in the latter 
scheme the system reached the steady state until 10 5 MC steps. In Fig. H] we plot the 
quantity p versus the time t (in MC steps) for all procedures after discarding 1 x 10 6 initial 
MC steps. We considered an initial configuration filled with particles and the densities are 
evaluated each 3 x 10 5 MC steps. With the exception of the arithmetic 1/T criterion, where 
the simulation gets trapped in the liquid phase (p ~ 1) the whole time of simulation, the 
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FIG. 4. The density p versus the time t at the phase coexistence (D* , 7\)=(8. 0000, 0.5000) after 
discarding 10 6 initial MC steps for the CEM, ad hoc, geometric and arithmetic 1/T distributions 
(graphs (a), (6), (c) and (d), respectively). The tie lines in graphs (a), (b), (c) denotes the stationary 
Po- 



system is able to cross the free-energy barriers properly in the other cases. This can be 
viewed by the fluctuations around its equilibrium value po, whose density averages p are 
consistent with po = 2/3 for arithmetic T, geometric and CEM. In the next application, 
we shall see that the choice of the temperature interval will have more influence on the 
tunneling. 



C. Bell-Lavis model 

In the last part of this paper, we study the BL model in triangular lattice of size L = N = 
18 using periodic boundary conditions. First, we repeat the analysis performed for the Ising 
model in the LDL-HDL second-order transition. We recall that the density of particles p is 
not the order-parameter 0, since p ^ in both liquid phases. A previous study [24] showed 
that the appropriate <fi is the difference between the fullest pi and the emptiest pj density 
sublattices given by (ft = pi — pj. In Fig. [5] we plot the auto-correlation functions and C u 
for all distributions. In particular, for the critical point located at (p, c , T c ) = (— 1.000, 0.430), 
we distributed R — 2 = 4 replicas between T\ = T c and T R = 0.730, whose entropies per 
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CEM ad hoc 

0.4300 0.430 

0.4601 0.468 

0.5022 0.518 

0.5589 0.576 

0.6340 0.645 

0.7300 0.730 

TABLE III. Temperature set of R = 6 replicas for the BL model at £t c = —1.000 obtained for the 
CEM and ad hoc distributions. 

0.8 
0.6 
0.4 
0.2 

°0 20 40 60 80 100 "0 20 40 60 80 100 
X X 

FIG. 5. Time displaced auto-correlation functions and C u versus time r (in MC steps) for the 
BL model at the critical point (fi c ,T c ) = (—1.000,0.430) for N = L = 18 and different temperature 
schedules. Results for the arithmetic 1/T have been omitted as they are similar to the ad hoc ones. 

site are si/ks = 0.5550(2) and sn/ks = 0.8763(2), respectively. We have also considered 
exchanges between every adjacent and non-adjacent (here between every second and third) 
temperatures. Table ITTT1 shows the temperature set for the CEM and ad hoc cases. As in the 
Ising model, and C u also decay faster at the critical point when temperatures are chosen 
using the CEM. Repetition for other critical points leads to the same conclusion. 

In addition to the previous study, we also investigate the first-order phase transition gas- 
LDL occurring at low temperatures. Numerical simulations have been carried out at T\ = 
0.1000. For this temperature, the phase transition takes place at ft* = —1.6500(1), which 
is identical (up to the fourth decimal level) to the transition point ft = —1.65 calculated 
at T = 0. In Fig. [61 we plot the time evolution of the density of molecules p = (crj) 
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FIG. 6. Time decay of the order parameter p versus time (in MC steps) for the BL model at 
the phase coexistence (p* , T\)=(— 1.6500, 0.1000) for N = L = 18 and different distribution of 
temperatures. The tie line for p = 1/2 denotes its stationary value pq. In the inset we plot, for the 
ad hoc and CEM schedules, the decay of p considering only adjacent swaps. 



starting from an initial configuration filled by molecules. We consider extreme temperatures 
T\ = 0.1000 and Tr = 0.4200, with corresponding entropies per site given by s\jks = 1CT 5 
and sn/kB = 0.4604(1), respectively. By considering in all cases a set of R = 6 replicas 
we have, for the CEM case, the intermediate temperatures T2 = 0.2837, T3 = 0.3456, 
T 4 = 0.3756 and T 5 = 0.3973. As for the BEG model, with the CEM the system crosses 
the entropic barriers more frequently than with other criteria, which can be identified by 
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29|. On the other 



the faster convergence of p toward its equilibrium value po ~ 1/2 
hand, for the other procedures, the system remains a larger number of MC steps trapped in 
metastable configurations and until 3 x 10 5 MC steps the density has not yet converged to po- 
When we consider only adjacent replica exchanges, the system gets trapped in metastable 
states for all distributions. In the inset of Fig. |6] we plot the decay of p only for the CEM 
and ad hoc schedules. 

We also show in Fig. [7] the density p versus t also starting from an initial configuration 
filled by molecules after discarding 10 6 initial MC steps. In the BL model, only with the 
CEM and ad hoc criteria the system crosses frequently the free-energy barriers, although the 
tunneling is rather more frequent with the CEM than the ad hoc distribution. Only in these 
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0.8 - 

. 6 1 — , — i — , — i — , — i — , — i — , — i — , — : 

50 100 150 200 250 300 

t(3.10 5 MC steps) 

FIG. 7. The density p versus the time t at the phase coexistence (p,*, 2\)=(— 1.6500, 0.1000) after 
discarding 10 6 initial MC steps for the CEM, ad hoc, geometric and arithmetic 1/T distributions 
(graphs (a),(6),(c) and (d), respectively) and N = L = 18. The tie lines in graphs (a),(6),(c) 
denotes the stationary pq. 

cases the tunneling between the different phases give an equilibrium value p = 0.504(3), 
consistent with po — 1/2- On the other hand, for arithmetic T (not shown) and geometric 
schedules the tunneling is much less frequent than the CEM ones, whereas for the arithmetic 
1/T the systems gets trapped the whole simulation in a metastable state generated by the 
initial configuration. 

As a final remark, it is worth mentioning that the observed differences in the results 
yielded by the various temperature schemes are less pronounced in discontinuous transitions 
occurring at high temperatures. In addition, by increasing both the number of replicas and 
the degree of non-adjacent exchanges the results also become closer. 

V. DISCUSSION AND CONCLUSION 

We presented a comparative study among different protocols for the choice of the tem- 
perature set in the parallel tempering method. We focused our attention on five criteria 
denoted arithmetic, geometric progressions in the temperature, arithmetic progression in 
the inverse temperature, ad hoc distribution and the constant entropy method (CEM). In 
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this last case, we considered an alternative direct MC method for evaluating the entropy 
which avoids numerical integrations of the specific heat or other thermodynamic quantities. 
We have considered rather few number of replicas (R — 6 for continuous and discontinu- 
ous phase transitions) and adjacent and non- adjacent replica exchanges. Different systems 
undergoing first and second-order phase transitions have been undertaken. In all cases, 
the temperature selection via the difference of entropy method revealed more advantageous. 
More specifically, at the criticality (where configurations generated by standard algorithms 
become strongly correlated) the time displaced correlation functions decay faster when tem- 
peratures are chosen with the CEM. This behavior can be understood that near criticality 
(where a small change of temperature provoke a large change of entropy) the CEM gives 
more concentrated intermediate temperatures than all distributions. Thus, replicas at the 
lowest temperature (T — 7\) display a larger probability of exchanging configurations than 
the other cases. Since the time correlation decays faster for T > T c than T = T c , the more 
frequent exchanges provide the system at Ti decays faster. For discontinuous transitions at 
low temperatures, where high entropic barriers do not allow the system to cross the phase 
frontiers properly (also when simulated by conventional algorithms) and hence the choice of 
the adjacent temperatures may play a crucial role, the CEM has also offered a rather effi- 
cient recipe for determining the temperature set. Within the CEM the lower temperatures 
are more sparse than with other schemes and, though unlikely, a successful replica exchange 
allows the system to evolve to configurations which are able to cross the high free energy 
barriers faster than other distributions. We have also distributed temperatures following 
an ad hoc scheme, in such a way that the exchange probability between adjacent replicas 
was about 30%. Although this method has shown to be more efficient than arithmetic and 
geometric schedules at the phase transition, it is inferior than the CEM. In summary, our 
comparative study effects the CEM as an useful tool for obtaining the temperature sched- 
ule to be used in numerical simulations of phase transitions through the parallel tempering 
method. 
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